package edu.spbsu.nonsm00th.qdiff;

import edu.spbsu.nonsm00th.computgeom.chull.ConvexHull;
import edu.spbsu.nonsm00th.computgeom.util.ConvexHullUtil;
import edu.spbsu.nonsm00th.model.Point;
import edu.spbsu.nonsm00th.qdiff.parser.RPNTranslator;
import edu.spbsu.nonsm00th.qdiff.parser.exception.ParserException;
import edu.spbsu.nonsm00th.qdiff.parser.word.Argument;
import edu.spbsu.nonsm00th.qdiff.parser.word.Function;
import edu.spbsu.nonsm00th.qdiff.parser.word.Unit;
import edu.spbsu.nonsm00th.qdiff.parser.word.Word;
import edu.spbsu.nonsm00th.qdiff.parser.word.Number;
import edu.spbsu.nonsm00th.qdiff.parser.word.Operation;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;
import java.util.Stack;

/**
 * Quasidifferential calculator class.
 * 12.10.2011 - Only for dimension 2 for now.
 * 
 * @author mbuzdov
 */
public class Calculator {
    private Stack<Quasidiff> qdiff = new Stack<Quasidiff>();
    private Stack<Double> val = new Stack<Double>();
    private int dimension;

    public Quasidiff calculate(String expression, Point point) throws ParserException {
        dimension = point.getDimension();
        if (2 != dimension) {//TODO: implement for other dimensions - just remove this.
            throw new IllegalStateException("Quasidiff calculator implemented only for 2nd dimension for now.");
        }
        RPNTranslator translator = new RPNTranslator();
        List<Unit> rpn = translator.translate(expression);
        verifyDimension(rpn);
        for (Unit next : rpn) {
            Word word = next.getWord();
            if (next.isNumber()) {
                List<Point> subdiff = new ArrayList<Point>(Arrays.asList(new Point(dimension)));
                List<Point> superdiff = new ArrayList<Point>(Arrays.asList(new Point(dimension)));
                qdiff.push(new Quasidiff(subdiff, superdiff));
                val.push(((Number)word).getValue());
            } else if (next.isArgument()) {
                int i = ((Argument)word).getI();
                Point p = new Point(dimension);
                p.setXi(i, 1.0);
                List<Point> subdiff = new ArrayList<Point>(Arrays.asList(p));
                List<Point> superdiff = new ArrayList<Point>(Arrays.asList(new Point(dimension)));
                qdiff.push(new Quasidiff(subdiff, superdiff));
                double v = point.getXi(i);
                val.push(v);
            } else {
                try {
                    if (next.isOperation()) {
                        process((Operation)next.getWord(), next.getArgumentsCount());
                    } else {//next.isFunction()
                        process((Function)next.getWord(), next.getArgumentsCount());
                    }
                } catch(Exception e) {
                    //TODO: add exception customizations and remove this try/catch block
                    throw new IllegalArgumentException("There is division by zero or another inappropriate operation.");
                }
            }
        }
        Quasidiff result = qdiff.pop();
        //TODO: find convex hulls by appropriate class and in appropriate place!
        ConvexHull chull = ConvexHullUtil.getD2ConvexHull();
        chull.convexHull(result.getSubdiff());
        chull.convexHull(result.getSuperdiff());
        return new Quasidiff(result.getSubdiff(), result.getSuperdiff());
    }

    private void process(Operation operation, int argCount) {
        if (Operation.ADDITION == operation) {
            if (2 == argCount) {
                qdiff.push(qdiff.pop().add(qdiff.pop()));
                val.push(val.pop() + val.pop());
            }
        } else if (Operation.SUBTRACTION == operation) {
            if (1 == argCount) {
                qdiff.push(qdiff.pop().multiply(-1.0));
                val.push(-val.pop());
            } else {//2 == argCount
                qdiff.push(qdiff.pop().multiply(-1.0).add(qdiff.pop()));
                val.push(-val.pop() + val.pop());
            }
        } else if (Operation.MULTIPLICATION == operation) {
            Quasidiff q2 = qdiff.pop();
            Quasidiff q1 = qdiff.pop();
            Double f2 = val.pop();
            Double f1 = val.pop();
            qdiff.push(q1.multiply(f2).add(q2.multiply(f1)));
            val.push(f1 * f2);
        } else if (Operation.DIVISION == operation) {
            Double f2 = val.pop();
            if (0 == f2) {
                throw new IllegalStateException("Dividing by zero!");//TODO: customize exceptions
            }
            Double f1 = val.pop();
            Quasidiff q2 = qdiff.pop();
            Quasidiff q1 = qdiff.pop();
            qdiff.push(q1.multiply(1.0/f2).add(q2.multiply(-f1/(f2*f2))));
            val.push(f1/f2);
        } else {//Operation.EXPONENTIATION == operation
            /* f(x)^g(x) = e^{g*ln(f)} */
            double g = val.pop();
            double f = val.pop();
            Quasidiff Dg = qdiff.pop();
            Quasidiff Df = qdiff.pop();
            Quasidiff Dlnf = Df.multiply(1.0/f);
            Quasidiff Dglnf = Dg.multiply(Math.log(f)).add(Dlnf.multiply(g));
            qdiff.push(Dglnf.multiply(Math.pow(f, g)));
            val.push(Math.pow(f, g));//TODO: pow rules?
        }
    }

    private void process(Function function, int argCount) {
        //TODO: ? incapsulate differential operation and replace this code with D(f(g(x)))=df*D(g(x))
        //TODO: maybe better to refactor with switch
        //TODO: verify function continuous
        //TODO: ? use Math.toRadians(f)
        if (Function.SIN == function) {
            double f = val.pop();
            qdiff.push(qdiff.pop().multiply(Math.cos(f)));
            val.push(Math.sin(f));
        } else if (Function.COS == function) {
            double f = val.pop();
            qdiff.push(qdiff.pop().multiply(-Math.sin(f)));
            val.push(Math.cos(f));
        } else if (Function.TG == function) {
            double f = val.pop();
            qdiff.push(qdiff.pop().multiply(1.0/(Math.cos(f)*Math.cos(f))));
            val.push(Math.tan(f));
        } else if (Function.CTG == function) {
            double f = val.pop();
            qdiff.push(qdiff.pop().multiply(-1.0/(Math.sin(f)*Math.sin(f))));
            val.push(1.0/Math.tan(f));
        } else if (Function.ARCSIN == function) {
            double f = val.pop();
            qdiff.push(qdiff.pop().multiply(1.0/Math.sqrt(1.0 - f*f)));
            val.push(Math.asin(f));
        } else if (Function.ARCCOS == function) {
            double f = val.pop();
            qdiff.push(qdiff.pop().multiply(-1.0/Math.sqrt(1.0 - f*f)));
            val.push(Math.acos(f));
        } else if (Function.ARCTG == function) {
            double f = val.pop();
            qdiff.push(qdiff.pop().multiply(1.0/(1.0 + f*f)));
            val.push(Math.atan(f));
        } else if (Function.ARCCTG == function) {
            double f = val.pop();
            qdiff.push(qdiff.pop().multiply(-1.0/(1.0 + f*f)));
            val.push(Math.PI/2 - Math.atan(f));
        } else if (Function.EXP == function) {
            double f = val.pop();
            qdiff.push(qdiff.pop().multiply(Math.exp(f)));
            val.push(Math.exp(f));
        } else if (Function.LN == function) {
            double f = val.pop();
            qdiff.push(qdiff.pop().multiply(1.0/f));
            val.push(Math.log(f));
        } else if (Function.POWER == function) {
            /* pow{f(x),g(x)} = f(x)^g(x) = e^{g*ln(f)} */
            double g = val.pop();
            double f = val.pop();
            Quasidiff Dg = qdiff.pop();
            Quasidiff Df = qdiff.pop();
            Quasidiff Dlnf = Df.multiply(1.0/f);
            Quasidiff Dglnf = Dg.multiply(Math.log(f)).add(Dlnf.multiply(g));
            qdiff.push(Dglnf.multiply(Math.pow(f, g)));
            val.push(Math.pow(f, g));//TODO: pow rules?
        } else if (Function.LOG == function) {
            /* log{g(x),f(x)} = ln{f(x)} / ln{g(x)} */
            double f = val.pop();
            double g = val.pop();
            Quasidiff Df = qdiff.pop();
            Quasidiff Dg = qdiff.pop();
            Quasidiff Dlnf = Df.multiply(1.0/f);
            Quasidiff Dlng = Dg.multiply(1.0/g);
            qdiff.push(Dlnf.multiply(1.0/Math.log(g)).add(Dlng.multiply(-Math.log(f)/(Math.log(g)*Math.log(g)))));
            val.push(Math.log(f)/Math.log(g));
        } else if (Function.SH == function) {
            double f = val.pop();
            qdiff.push(qdiff.pop().multiply(Math.cosh(f)));
            val.push(Math.sinh(f));
        } else if (Function.CH == function) {
            double f = val.pop();
            qdiff.push(qdiff.pop().multiply(Math.sinh(f)));
            val.push(Math.cosh(f));
        } else if (Function.TH == function) {
            double f = val.pop();
            qdiff.push(qdiff.pop().multiply(1.0/(Math.cosh(f)*Math.cosh(f))));
            val.push(Math.tanh(f));
        } else if (Function.CTH == function) {
            double f = val.pop();
            qdiff.push(qdiff.pop().multiply(-1.0/(Math.sinh(f)*Math.sinh(f))));
            val.push(1.0/Math.tanh(f));
        } else if (Function.MAX == function) {
            /* max{fi(x)}, i=1,..,N} */
            double[] fi = getArgs(argCount);
            List<Integer> actives = getMaxActive(fi);
            val.push(fi[actives.get(0)]);
            qdiff.push(getMax(getQuasidiffs(argCount), actives));
        } else if (Function.MIN == function) {
            /* min{fi(x)}, i=1,..,N} */
            double[] fi = getArgs(argCount);
            List<Integer> actives = getMinActive(fi);
            val.push(fi[actives.get(0)]);
            qdiff.push(getMin(getQuasidiffs(argCount), actives));
        }
    }

    /* TODO: incapsulate/refactor this part, add SetUtil class */
    private Point sum(Point p1, Point p2) {
        double[] x = new double[p1.getDimension()];
        for (int i = 0; i < p1.getDimension(); i++) {
            x[i] = p1.getXi(i) + p2.getXi(i);
        }
        return new Point(x);
    }

    private List<Point> sum(List<Point> l1, List<Point> l2) {
        List<Point> sum = new ArrayList<Point>();
        for (Point p1 : l1) {
            for (Point p2 : l2) {
                sum.add(sum(p1, p2));
            }
        }
        return sum;
    }

    private List<Point> multiply(List<Point> list, double c) {//TODO: refactor this
        List<Point> result = new ArrayList<Point>();
        for (Point p : list) {
            Point temp = new Point(p.getDimension());
            for (int i = 0; i < p.getDimension(); i++) {
                temp.setXi(i, p.getXi(i) * c);
            }
            result.add(temp);
        }
        return result;
    }

    private Quasidiff getMax(Quasidiff[] quasis, List<Integer> actives) {
        List<Point> subdiff = new ArrayList<Point>(Arrays.asList(new Point(2)));
        List<Point> superdiff = new ArrayList<Point>(Arrays.asList(new Point(2)));
        for (int next : actives) {
            superdiff = sum(superdiff, quasis[next].getSuperdiff());
            List<Point> npart = quasis[next].getSubdiff();
            for (int i : actives) {
                if (next != i) {
                    npart = sum(npart, multiply(quasis[i].getSuperdiff(), -1.0));
                }
            }
            subdiff.addAll(npart);
        }
        return new Quasidiff(subdiff, superdiff);
    }

    private Quasidiff getMin(Quasidiff[] quasis, List<Integer> actives) {
        List<Point> subdiff = new ArrayList<Point>(Arrays.asList(new Point(2)));
        List<Point> superdiff = new ArrayList<Point>(Arrays.asList(new Point(2)));
        for (int next : actives) {
            subdiff = sum(subdiff, quasis[next].getSubdiff());
            List<Point> npart = quasis[next].getSuperdiff();
            for (int i : actives) {
                if (next != i) {
                    npart = sum(npart, multiply(quasis[i].getSubdiff(), -1.0));
                }
            }
            superdiff.addAll(npart);
        }
        return new Quasidiff(subdiff, superdiff);
    }

    private double[] getArgs(int argCount) {
        double[] result = new double[argCount];
        for (int i = 0; i < argCount; i++) {
            result[argCount - 1 - i] = val.pop();
        }
        return result;
    }

    private Quasidiff[] getQuasidiffs(int argCount) {
        Quasidiff[] result = new Quasidiff[argCount];
        for (int i = 0; i < argCount; i++) {
            result[argCount - 1 - i] = qdiff.pop();
        }
        return result;
    }

    private List<Integer> getMaxActive(double[] values) {
        return getActive(values, true);
    }

    private List<Integer> getMinActive(double[] values) {
        return getActive(values, false);
    }

    private List<Integer> getActive(double[] values, boolean isMax) {
        List<Integer> result = new ArrayList<Integer>();
        double m = values[0];
        result.add(0);
        for (int i = 1; i < values.length; i++) {
            double next = values[i];
            if (next == m) {
                result.add(i);
            } else if (isMax && next > m || !isMax && next < m) {
                result.clear();
                result.add(i);
                m = values[i];
            }
        }
        return result;
    }

    private void verifyDimension(List<Unit> list) {
        for (Unit next : list) {
            if (next.isArgument() && ((Argument)next.getWord()).getI() >= dimension) {
//                int i = ((Argument)next.getWord()).getI();
                throw new IllegalStateException("There is argument in expression, which dimension is more than possible one.");
                //TODO: customize exception: CalculatorException or etc.
            }
        }
    }
}
